# plot sunspot deviation means
# http://solarscience.msfc.nasa.gov/SunspotCycle.shtml

sunspots <- function()
{
	X <- read.table('data/spot_num.txt', header=TRUE)
	X[,3]
}


sunspots_statistics <- function(X)
{
	wmean <- rep(0, 133)
	wsd <- rep(0, 133)
	for(t in 1:133) {
		Xt <- X[c(t+133*(0:floor(length(X))))]
		wmean[t] <- mean(Xt, na.rm=TRUE)
		wsd[t] <- sd(Xt, na.rm=TRUE)
	}
	data.frame(mean=wmean, sd=wsd)
}


#par(mfrow=c(1,2), mar=c(1.8,2,0.2,0))

#################### MEAN

sunspots_means <- function(X)
{
wmean <- sunspots_statistics(X)$mean

plot(1:133, wmean, pch=20, main=NULL, xlab='month', ylab='sunspot numbers', ylim=c(20,80))
axis(side=2, at=seq(20,80,20))
grid(NULL, NULL, lwd=1)

# regression
# <65  and  >73
# we could also see R^2 from the model object

t <- 1:64
d <- data.frame(x=t, y=wmean[t])
m <- lm('y~x', d)
abline(m$coefficients[1], m$coefficients[2], col='red', lwd=2)

t <- 74:133
d <- data.frame(x=t, y=wmean[t])
m <- lm('y~x', d)
abline(m$coefficients[1], m$coefficients[2], col='red', lwd=2)
}


#################### STANDARD DEVIATION

sunspots_deviations <- function(X)
{
wsd <- sunspots_statistics(X)$sd

plot(1:133, wsd, pch=20, main=NULL, xlab='month', ylab='sunspot numbers', ylim=c(20,60))
axis(side=2, at=seq(20,60,10))
grid(NULL, NULL, lwd=1)

# first harmonic model
# see http://stats.stackexchange.com/questions/60500/how-to-find-a-good-fit-for-semi-sinusoidal-model-in-r

t <- 1:133
d <- data.frame(x=t, y=wsd)
period <- 133
m <- lm(y ~ sin((2*pi/period)*x) + cos((2*pi/period)*x), d)
lines(t, m$fitted, col='red', lwd=2)
}


